# Gelman Carlin retrospetive power and beyond-power analysis:

# Make sure directory is set to be the correct location
# (e.g., use setwd() to set directory to the LDAR_RCT_replication_package folder)

library(dplyr)
library(readr)

# Reads in outputs created by from site_level_for_analysis.do
param_loc <- file.path('Data/Intermediate/power_gelman_carlin_inputs.csv')
param_df <- read_csv(param_loc)


# Original Gelman-Carlin function:
retrodesign <- function(A, s, alpha=.05, df=Inf, n.sims=10000){
  z <- qt(1-alpha/2, df)
  p.hi <- 1 - pt(z-A/s, df)
  p.lo <- pt(-z-A/s, df)
  power <- p.hi + p.lo
  typeS <- p.lo/power
  estimate <- A + s*rt(n.sims,df)
  significant <- abs(estimate) > s*z
  exaggeration <- mean(abs(estimate)[significant])/A
  return(list(power=power, typeS=typeS, exaggeration=exaggeration))
}

# Augmented for the case where we don't know standard deviations 
#   and where standard deviations may vary across treatment/control
retrodesign_2 <- function(A, sd1, sd2, n1, n2, alpha = 0.05, n.sims=10000) {
  # Standard deviation for case where we estimate from sample standard deviation
  s <- ((sd1^2)/n1 + (sd2^2)/n2)^0.5
  # Satterthwaite degrees of freedom
  df <- ((sd1^2)/n1 + (sd2^2)/n2)^2 / 
    ((((sd1^2)/n1)^2)/(n1-1) + (((sd2^2)/n2)^2)/(n2-1))
  
  z <- qt(1-alpha/2, df)
  p.hi <- 1 - pt(z-A/s, df)
  p.lo <- pt(-z-A/s, df)
  power <- p.hi + p.lo
  typeS <- p.lo/power
  estimate <- A + s*rt(n.sims,df)
  significant <- abs(estimate) > s*z
  exaggeration <- mean(abs(estimate)[significant])/A
  return(list(df=df, s=s, power=power, typeS=typeS, exaggeration=exaggeration))
  
}

perc_size <- c(1, 0.5, 0.25, 0.1)

# Creates blank matrix for output:
beyond_power_mat <- matrix("", nrow = 9, ncol = 10)

# Naming rows and columns
beyond_power_mat[2, 1] <- "\\hline"
beyond_power_mat[6, 1] <- "\\hline"
beyond_power_mat[1,2] <- "Est eff."
beyond_power_mat[1,3] <- "Hypo. eff."
beyond_power_mat[1,4] <- "Stdev.\\textsubscript{C}"
beyond_power_mat[1,5] <- "Stdev.\\textsubscript{T}"
beyond_power_mat[1,6] <- "N\\textsubscript{C}"
beyond_power_mat[1,7] <- "N\\textsubscript{T}"
beyond_power_mat[1,8] <- "Power"
beyond_power_mat[1,9] <- "Type S"
beyond_power_mat[1,10] <- "Exag"

# Now fills in the matrix:
k <- 2
for (i in c(1,2)) {

  # Input components:  
  est_effect_size <- param_df$muT[i] - param_df$muC[i]
  sdC <- param_df$sdC[i]
  sdT <- param_df$sdT[i]
  nC <- param_df$nC[i]
  nT <- param_df$nT[i]

  for (j in seq_along(perc_size)) {
    
    # scalar for effect size:
    perc_size_j <- perc_size[j]
    hypothesized_true_effect <- perc_size_j * est_effect_size
  
    # Runs the power calculation:    
    design_results <- retrodesign_2(hypothesized_true_effect, sdC, sdT, nC, nT)
    
    # Outputs results to matrix:
    beyond_power_mat[k, 2] <- round(est_effect_size, digits=2) %>% sprintf("%.2f", .)
    beyond_power_mat[k, 3] <- round(hypothesized_true_effect, digits=2) %>% sprintf("%.2f", .)
    beyond_power_mat[k, 4] <- round(sdC, digits=2) %>% sprintf("%.2f", .)
    beyond_power_mat[k, 5] <- round(sdT, digits=2) %>% sprintf("%.2f", .)
    beyond_power_mat[k, 6] <- nC
    beyond_power_mat[k, 7] <- nT
    beyond_power_mat[k, 8] <- round(design_results$power, digits=2) %>% sprintf("%.2f", .)
    beyond_power_mat[k, 9] <- round(design_results$typeS, digits=2) %>% sprintf("%.2f", .)
    beyond_power_mat[k, 10] <- round(design_results$exaggeration, digits = 2) %>% sprintf("%.2f", .)

    # Moves to next k:
    k <- k + 1
    
  }
}

beyond_power_mat 

# outputting to Latex
outputtable_to_LaTeX <- function(arr) {
  rows <- apply(arr, MARGIN = 1, paste, collapse = " & ")
  matrix_string <- paste(rows, collapse = " \\\\ \n ")
  return(paste("\\begin{tabular}{| l c c c c c c | c c c |}\\hline\\hline \n", 
               matrix_string, "\\\\ \n",
               "\\hline \\hline \\end{tabular}"))
}

# Exports:
beyond_power_mat_path <- file.path("Paper/Figures/gelman_carlin_beyond_power.tex")
cat( outputtable_to_LaTeX(beyond_power_mat), file = beyond_power_mat_path)


